Quantum rainbows and ergodicity in many-body dynamics 
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We consider rainbows in Fock space that form following a quench in a Bose Josephson junction. 
In the Gross-Pitaevskii mean-field theory the rainbows are singular caustics, but in the second- 
quantized theory they are described by well behaved Airy functions. We discuss the analyticity 
of the macroscopic limit and also show that the long-time dynamics are ergodic. Our results are 
relevant to the question posed by Berry [M.V. Berry, Nonlinearity 21, T19 (2008)]: are there 
circumstances when it is necessary to second-quantize wave theory in order to avoid singularities? 



In the geometrical theory of light, rainbows are caus- 
tics formed by the focusing of the Sun's rays by rain- 
drops. Caustics are the envelopes of families of rays and 
the light intensity diverges upon them: they are catas- 
trophes where the ray theory fails. In order to tame 
the divergences it is necessary to go to the next level 
of sophistication, namely the wave theory of light. For 
rainbows this problem was solved by Airy in 1838 PJ. 

In this paper we investigate caustics that form in Fock 
space during quantum many-body (MB) dynamics. The 
analog of ray and wave theory are mean-field (MF) and 
MB theory, respectively. Caustics in Fock space are 
places where a classical field must be second quantized 
in order to remove divergences. Just as wave theory in- 
troduces the notion of phase, MB theory introduces sec- 
ond quantization, i.e. discreteness of quanta, and this 
smooths the caustic singularity. There are many situa- 
tions where MF theory is quantitatively inaccurate, but 
at its caustics it fails qualitatively, predicting infinities, 
and second quantization is forced upon us. 

Related ideas have previously been discussed in the 
context of phase singularities in light waves [2 By 
contrast, our Fock space caustics are amplitude singular- 
ities. Caustics also occur in atom optics, e.g. during the 
dynamical diffraction of atoms by a standing wave of light 
[4] as observed by an experiment using Bose-Einstein con- 
densates (BECs) [5] , and during the expansion of a BEC 
[5]. However, these latter studies concerned caustics in 
real or momentum space where matter wave interference 
from first quantization smooths the caustic singularity. 

Cold atoms are useful for studying quantum MB dy- 
namics because of their long coherence times and high 
tunability |7]. A key experimental workhorse in such 
studies is the degenerate Bose gas in an optical lattice 
which can realize strongly interacting (non MF) 
states, such as Mott insulators [14]. Non-equilibrium 
states can be generated by suddenly changing the lattice 
depth [HHH] , and adiabatic MB dynamics have also been 
studied |13j . The system we consider here is a poor man's 
optical lattice, namely a Bose Josephson junction (BJJ) 
made of two lattice sites. Despite its relative simplicity, 
the BJJ can be used as a qualitative guide to the many- 
site problem [9]. Furthermore, its MB dynamics can be 
understood analytically. Two pioneering cold atom ex- 
periments have demonstrated Josephson effects in the 



BJJ: plasma oscillations and macroscopic self-trapping 
[T51 [TB] (a manifestation of the ac Josephson effect) were 
observed in [T7] , and both the ac and dc Josephson effects 
were observed in [18] . The height of the tunneling barrier 
can be controlled, and the effect this has on Josephson 
oscillations has been measured in [19j . 

There is a large body of theoretical work on the BJJ. 
Of particular relevance is [TS] , which pointed out the pos- 
sibility of a collapse and revival sequence in the distribu- 
tion of atom number differences between the two sites if 
all the atoms are initially localized on just one site. In 
this work we consider the collapses and revivals that fol- 
low a quench in the tunneling barrier and show that each 
revival corresponds to the birth of a MB rainbow. At 
long times Fock space becomes thick with rainbows. 

The Hamiltonian for a BJJ with TV bosons is [T6ll20l[2T] 

H = ^EcU^ - EJ^yl - 4n^ /N^ coscj) (1) 

where n — {Ni — Nr)/2 is half the number difference be- 
tween the left- and right-hand BECs, and (f) = (j)^ — ipi 
is the phase difference. In the MF description, which 
is provided by the Gross-Pitaevskii theory [20], n and (j) 
are numbers taking on continuous and well-defined val- 
ues (like the amplitude and phase of the electric field in 
Maxwell's theory of light) whereas in the MB description 
they are operators obeying [0, n] = i, which forces their 
eigenvalues to be discrete. For notational simplicity, in 
the latter case we use n interchangeably to represent both 
the operator and its eigenvalues. Cold atom experiments 
have direct access to n and (j) (but not at the same time) . 

The charging energy Ec derives from interatomic in- 
teractions and Ej from tunneling between the sites. The 
ratio A = 2Ej/Ec can be used to define three regimes 
[55]: Fock A < 1, Josephson 1 < A < A^^, and Rabi 
A ^ N"^. In the Josephson and Fock regimes (A <g; N^), 
and at low energies (i.e. when n <^ N), the a/1 — 4rfi /N"^ 
factor can be set to unity and the hamiltonian takes on 
the universal form of a rigid pendulum. This is the situ- 
ation we shall assume from now on. The phase space of 
the classical pendulum is divided by a separatrix (a con- 
tour with energy Ej corresponding to the barrier tops of 
the cos (j) potential) into two distinct dynamical regions: 
inside the separatrix the pendulum oscillates back and 
forth, and outside it the pendulum makes complete ro- 
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tations. The simplified hamiltonian still contains macro- 
scopic self trapping, which derives from rotational motion 
of the pendulum. Furthermore, due to the powerful con- 
cept of structural stability (see below), our results con- 
cerning rainbow caustics also apply qualitatively even in 
the Rabi regime. 

In the pendulum analogy the first term in the hamil- 
tonian is kinetic energy, and therefore oc Ti^ . Thus, 
A oc m [5S]. Small values of A correspond to the 
quantum regime where there are few energy eigenstates 
below the separatrix, and large values to the semiclassi- 
cal regime where there are many. In fact, in the Bose- 
Hubbard (BH) model X (x N (although there are other 
possibilities beyond the BH model [231 [H]), and we can 
identify 1/y/N as playing the role of Planck's constant in 
the MB problem. Because we are interested in how sin- 
gularities develop as N becomes large, our calculations 
will be in the semiclassical regime. 

Prior to the quench we assume that the tunneling bar- 
rier is high, and the state of the B JJ can be approximated 
by the Fock state n = 0. This is the ground state when 
A = 0. However, we can once again invoke the structural 
stability of rainbow caustics to assert that, even if the 
initial state corresponds to a small distribution of Fock 
states, our results will remain qualitatively correct. At 
t — the tunneling barrier is suddenly lowered so that 
the BJJ goes from the Fock to the Josephson regime. 

Consider the MF dynamics first. While each MF so- 
lution has a perfectly defined phase, the initial MB state 
corresponds to a superposition of all possible phases, and 
is highly non-classical. To mimic this we can form the 
sum of a large number of "rays" , a selection of which are 
shown in the lower part of Fig. [T] with initial phases (f)o 
uniformly distributed over [0, 27r). The MF rays obey 

n{t) = y2Asin[0o/2]cn{a;pii-f K(sin^[0o/2])|sin^[0o/2]} 

where Wpi = y/EcEj/h is the plasma frequency which 
gives the oscillation frequency of low-lying excitations, 
and cn{u|u} and K{v) are a Jacobi elliptic function and 
a complete elliptic integral of the first kind, respectively, 
|28j . The cosine potential in the Hamiltonian is not a 
perfect lens due to its anharmonicity, but instead gives 
rise to extended caustics that proliferate with time. Each 
new caustic is born as a cusp point (three rays focused to 
the same point) at the times tm — rmr/ujpi, where m = 
0, 1, 2, . . ., and then evolves into a pair of rainbows (fold 
caustics due to the focusing of two rays) that gradually 
move out to n = ±\/2A. 

Classical dynamics is a gradient map (Fermat's prin- 
ciple), and the singularities of gradient maps are clas- 
sified by catastrophe theory [26l [27] . The classification 
is in terms of co-dimension k, = dimension of space — 
dimension of caustic. The only structurally stable caus- 
tics in two dimensions are cusp points (k — 2) and fold 
lines {k = 1), and these are exactly what we see in Fig. 
[T] which is a two-dimensional space formed by combining 
one-dimensional Fock space with time. 




FIG. 1. MB and MF dynamics of the BJJ following the 
quench. Time is measured in units of the inverse plasma fre- 
quency, and n is the Fock space coordinate (half the number 
difference of bosons between the two sites) . Above: MB prob- 
ability distribution (for A — 200). Below: MF rays, where 
each ray has a different initial phase difference (po, lying in 
the range — tt . . . tt in steps of 7r/45. 

Turning now to the MB dynamics, we express 
the wave function as an energy eigensum |^(t)) = 
J2j o^jlEj) eyipl—iEjt/h]. If we write the energy eigen- 
states as sums over Fock states \Ej) — X^n^j.nl'^)' 
then, by the sudden approximation, the coefficients in 
the eigensum following the quench are aj = Aj^q because 
the initial state is n = 0. In fact, the energy eigenstates of 
the rigid quantum pendulum are Mathieu functions |28j . 
In the semiclassical regime, the eigenstates above the sep- 
aratrix can be dropped from the sum because these states 
are rotor states corresponding to finite values of "angu- 
lar momentum" n, and have only an exponentially small 
amplitude at n = 0, meaning that aj ~ 0. 

Buried in the eigensum representation for |^(i)) is the 
characteristic Airy function structure of the two rain- 
bows that are visible in Fig. [2| This figure gives a slice 
through Fig. [T] corresponding to the red line. In order 
to bring out the rainbows explicitly we follow Berry ^\ 
and apply the Poisson resummation formula X]j=o /(■?) ~ 
Em=-oo ir exp[2'Kimj]dj , where /(j) refers to the 
j'^ term in the energy eigensum. This transformation, 
which is exact, converts the sum over eigenstates into 
a path integral representation which sums over families 
of classical paths (rays) 'SOJ. Each family, labelled by 
the Maslov index m, has geometric significance and cor- 
responds to a rainbow pair. Referring to Fig. [T] we see 
that the MF caustics correspond to periodic collapses and 
revivals in the MB wave function. 

The number of states M we need to include in the 
eigensum is the number below the separatrix, and in the 
semiclassical regime Af w •\/2A [31. In the BH model 
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oc y/N. While this is a small number in comparison to 
total number of states + 1, it still diverges as N ^ co. 
By contrast, only two terms are needed in the Poisson 
sum for the time slice shown in Fig. |2] 

We now give a sketch of the Poisson resummation of 
the energy eigensum (the details will be published else- 
where). Firstly, we replace the Mathieu functions \Ej) 
by their WKB approximations We then apply the 
Poisson summation formula to obtain 



1 

2^ 



E 



/AC 



d/3 

(3) 

where the denominator D ~ (1 — (j/^ — — /?^)^''*, 

and we have introduced y — n/^/X and /3 — E/\, which 
are classical versions of the number difference and energy 
eigenvalues. The energies below the separatrix lie in the 
range — 1 < /3 < 1. The four phases are given by 
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where Pi{y,l3) = 
2VrT73E{arccos[?;2 



P2{y.P) 

and P3(/3) 



?/arccos[?/^ 
/3]/2|2/(l 

P2(0,/3). In these expressions E{u|u} is the Jacobi el- 
liptic integral of the second kind [28^ . 

The integrals appearing in Eq. ^ can be evaluated by 
the uniform approximation ^29) , which is a version of the 
method of stationary phase. When there is only a single 
stationary point /3i, the phase v, which represents either 
A, B, C, or I?, can be mapped onto a quadratic function 
which gives a Gaussian integral and the result 
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exp[i(VAz^(/3i)±V4)] 



(8) 



where the ± sign is fixed by the sign of j^"(/3i), which 
is the second derivative of the phase with respect to (3 
evaluated at the stationary point. The second derivative 
of phase A is A" = +Qi - (92 + Qa + (1 - 4m)(54, and 

in terms of this template the others are B" = -i 1 h, 

C" = — I 1 — , V" = — ! 1~. In these expressions 

Qi = F{g\h}/{2V2il + l3)}, Q2 = E{.g|/i}/{y2(l - /J^)}, 
Q3 = 2/(2/3 - y2)/{2(l - /?2)Vl-(y2_;3)2}, and Q, = 
[K{h) - 2E{tt / 2\h} /{I- 13)] /{2V2{1 + 13)}, where F{u\v} 
is the Jacobi elliptic integral of the first kind [28], and 
we have denoted g = arcsin[^ {I + (3 — y^)/(l + 13)] and 
h={l + f3)/2. 

On the bright side of a rainbow there are two station- 
ary points f3i and l32 which coalesce and annihilate on 
the fold line itself; on the dark side they become com- 
plex. Thus, on its bright side ^'rainbow displays two-wave 
interference fringes, and on the dark side it decays expo- 
nentially. Near a cusp there are three stationary points 




FIG. 2. Fock space rainbows 
distribution at the time t — 



number difference probabilit 



lity 
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37r/2aj~[ (red line in Fig. 
Due to the ±n symmetry of the solution, only the -|-ve haff 
shown. Dashed black curve: MF caustics given by summing 6 
million rays [Eq. ([2|]. Red dots: exact MB solution given by 
numerical diagonalization with A — 12000. Solid blue curve: 
semiclassical MB solution [Eqns ^ and with A = 12000. 
The semiclassical solution is only evaluated at integer values 
of n, but the points have been joined to guide the eye. 



which coalesce at the cusp point. For simplicity, we shall 
concentrate on rainbows (also, the cusps are confined to 
the immediate vicinity of n = 0). 

When there are two stationary points the phase can be 
mapped onto a cubic polynomial and the integral gives an 
Airy function Ai(C) = exp[i{s^/3 + Cs)]ds. Thus, 
the uniform approximation for ^E* due to a rainbow is 
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(9) 



where D = [v{pi) + v{p2)]/'2. and Av = 1/(^2) - viPi). 
The appearance of the derivative term Ai'(C) — dAi/dC 
in Eq. ( 9 1 is due to a series expansion of the amplitude 
of Eq. |3{ 29J. On the fold line (C = 0) the Airy term 
dominates the derivative term, but the latter makes a sig- 
nificant contribution away from the caustic. Note that 
the v" factors in the denominators precisely cancel di- 
vergences that would otherwise occur due to the use of 
WKB approximations for the Mathieu functions. 

As can be seen in Fig. [2] the semiclassical approx- 
imation gives an excellent match to the exact wave 
function (obtained by numerical diagonalization) when 
A = 12, 000. The only (real) stationary points we find 
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FIG. 3. Ergodicity in the BJJ at long times. Solid red curve: 
MF result for the number difference probability distribution, 
as given by Eq. (10 1, which is obtained by averaging over 
classical phase space .11. Blue dots: time average of the MB 
distribution over range Upit = 61.27r-62.27r with A = 12000. 
Dashed black curve: MB distribution at thermal equilibrium. 



at time t — (37r/2)ajp[ are in the m — and m — 1 
families, and so these are the only nonzero terms in the 
Poisson sum. The m = and m = 1 terms give the outer 
and inner rainbows, whose caustic singularities occur at 
ric/V^X « 0.97 and Uc/V^X « 0.52, respectively. The 
inner rainbow is not a clean Airy function because of in- 
terference with the m — term. However, for large A the 
Airy function near the caustic will dominate. From Eqns 
Q and ^ we find that, relative to its background, j'l'p 
evaluated on the fold caustic diverges as A^^^, see Table 
[ij At a cusp the relative divergence is X^^^ [17]. If we 
assume the BH model, the divergence in the macroscopic 
regime goes as N^^^ and iV^^"* for the folds and cusps, 
respectively. This shows that the macroscopic limit is 
non-analytic: in different regions of Fock space ^' has a 
different dependence on N. 





n <^ ric (bright) 


n — ric 


ric (dark) 




0{\-i) X cos[C>(A3)] 


0{\~i) 


0(A-3)exp[-C(A5)] 



TABLE I. Dependence of the number difference probability 
|3'(n)p on A as n is varied about a fold caustic (rainbow). 

Finally, let us consider the long time dynamics. In the 
MF theory, the ergodic average of the number difference 
probability distribution can be expressed as |31j : 



P{n) 



1 

2n 



dm 



K [m] ^y m(l — m){m — ri^)(l + rt-^ — m) 

(10) 

which is plotted in Fig. [3] along with the MB distribution, 
which has been averaged over one plasma oscillation pe- 
riod between the times t = (61.2-62.2)7rti;~[^, where there 
are 62 rainbows. Their close correspondence indicates 
that the MB dynamics is ergodic. This is to be expected 
even in this integrable system: there is only one degree 
of freedom and so the energy contour coincides with the 
phase space trajectory. For comparison, we have also in- 
cluded the MB distribution at thermal equilibrium (with 
the same average energy as the coherent state). The er- 
godic and thermal states are clearly distinguishable. 

Conclusions We have shown that Airy functions dec- 
orate the caustics that proliferate in Fock space follow- 
ing a quench in a BJJ. These second-quantized rainbows 
dominate the wave function in the macroscopic limit 
(N — >■ co) corresponding to the quantum-to-classical 
transition, which is non-analytic. These results are ro- 
bust against perturbations to the Hamiltonian and initial 
conditions because rainbows are structurally stable. The 
MB dynamics appear to be ergodic and agree with the 
MF result. It would be interesting to consider a greater 
number of lattice sites: Fock space then has more dimen- 
sions and higher order MB catastrophes are possible. 
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